A personalized and non-parametric framework for detecting changes in gait cycles
Saeb Ragani Lamooki1
Department of Mechanical & Aerospace Engineering, University at BuffaloJiyeon Kang2
Department of Mechanical & Aerospace Engineering, University at BuffaloLora A. Cavuoto3
Department of Industrial and Systems Engineering, University at BuffaloFadel M. Megahed4
Farmer School of Business, Miami UniversityAllison Jones Farmer5
Farmer School of Business, Miami UniversityJune 18, 2020
1 R Setup and Required Packages
In this project, the open-source R programming language is used to model/monitor changes in the gait cycles over time. R is maintained by an international team of developers who make the language available at The Comprehensive R Archive Network. Readers interested in reusing our code and reproducing our results should have R installed locally on their machines. R can be installed on a number of different operating systems (see Windows, Mac, and Linux for the installation instructions for these systems). We also recommend using the RStudio interface for R. The reader can download RStudio for free by following the instructions at the link. For non-R users, we recommend the Hands-on Programming with R for a brief overview of the software’s functionality. Hereafter, we assume that the reader has an introductory understanding of the R programming language.
In the code chunk below, we load the packages used to support our analysis.
if(require(checkpoint)==FALSE) install.packages("checkpoint") # check to see if checkpoint is installed; if not, install it
library(checkpoint) # package used to facilitate the reproducibility of our work
# a checkpoint of R packages on CRAN on June 08, 2020 to enable the reproduction of our work in the future
checkpoint("2020-06-08")
# check if packages are not installed; if yes, install missing packages
packages = c("tidyverse", "magrittr", # typical data analysis packages
"foreach", "parallel", # packages used for parallelizing some code chunks
"MALDIquant", # match closest points between two vectors
"MASS", # calculates the generalized inverse matrix
"R.matlab")
newPackages = packages[!(packages %in% installed.packages()[,"Package"])]
if(length(newPackages) > 0) install.packages(newPackages)
# using the library command to load all packages; invisible used to avoid printing all packages and dependencies used
invisible(lapply(packages, library, character.only = TRUE))
source("./Functions.R") # our custom built functions
set.seed(2020)
startTime <- Sys.time()2 Extracting, Transforming and Loading (ETL) the Gait Data into R
In this paper, we perform a secondary data analysis of the experimental data reported in Baghdadi et al. (2019), where fifteen participants were equipped with an IMU sensor wrapped on their right ankle. We used the IMU data to transform the linear acceleration signals from the local frame of reference to the global frame of reference. Next, we calculated the acceleration magnitude signal and smoothed it by applying a low pass butterworth filter with a cut-off frequency of 10 Hz. Further, we segmented the idividual gait cycles from foot-flat to foot-flat for each participant by tuning the parameters of the segmentation algorithm for each participant.
We load four different sets of data into R:
- Subject heights, which is read directly from the GitHub Repository of the Baghdadi et al. (2019) paper. Note that we have excluded Subject 13’s data due to ———–we need to insert some reasons here—————.
- Features, where we extract the x and y positions as well as stride duration and experimental time. The first three elements are used to compute stride length, height and duration which are three important kinematic features that can be used to model/monitor changes in gait. Note that we did not use this data directly from the Baghdadi et al. (2019) paper since our segmentation here was from foot-flat to foot-flat and a different frequency for the butterworth filter was used (given that the goals of our paper and the original paper are different).
- Profiles, where the points forming each gait segment are mantained. The profiles present an alternative approach to monitoring and modeling changes in the gait cycles.
- Video tags, where have coded when subjects enter and leave the main laboratory room. This information is used to ensure the consistency of the task being performed by the participant, i.e., we do not want the detected changes to correspond to a task change, but rather to indicate fatigue (see the discussion in our manuscript for more details).
The code chunk below provides our code for reading the three different datasets. The subject heights dataset is saved into a vector, which we use to scale the computed gait/stride length and height. The main outputs from this chunk are two lists (a) subjectFeaturesList which contains the scaled stride length, scaled stride height, stride duration and experimental time, and (b) subjectProfilesList which contains the sample IMU acceleration points making each acceleration profile.
# Reading the subject Heights from the Baghdadi et al. (2019) GitHub Repo:
# ------------------------------------------------------------------------
subject_heights <- read.csv("https://raw.githubusercontent.com/fmegahed/fatigue-changepoint/master/Data/Raw/DemographicData.csv") %>% pull(Height) %>% `*` (0.01) %>% .[-13] # pull height and convert into meters
# Computing gait features based on experimental data from Baghdadi et al. (2019):
# -------------------------------------------------------------------------------
featureFiles <- list.files(path = "../Data/matFiles/Features/", pattern = "*.mat", full.names = TRUE)
rawSubjectsData <- lapply(featureFiles, readMat) %>% # read all mat files in Feature List
lapply("[[", "M.i.k") # extract sublist titled "M.i.k" from all the lists
names(rawSubjectsData) <- str_extract_all(featureFiles, "Subject[:digit:]*", simplify = TRUE)
numRows <- lengths(rawSubjectsData)/17 # number of observations per Subject
subjectFeaturesList <- vector(mode = "list", length = length(rawSubjectsData)) # initialize list
names(subjectFeaturesList) <- names(rawSubjectsData)
for (i in seq(length(rawSubjectsData))) {
gaitLength <- rawSubjectsData[[i]][ 1:numRows[i] ] %>%
lapply(range) %>% lapply(diff) %>% unlist() %>%
`/` (subject_heights[i])
gaitHeight <- rawSubjectsData[[i]][ (1+numRows[i]):(2*numRows[i]) ] %>%
lapply(range) %>% lapply(diff) %>% unlist() %>%
`/` (subject_heights[i])
gaitDuration <- rawSubjectsData[[i]][ (1+ 2*numRows[i]):(3*numRows[i]) ] %>%
lapply(range) %>% lapply(diff) %>% unlist(recursive = TRUE)
expTime <- rawSubjectsData[[i]][ (1+ 16*numRows[i]):(17*numRows[i]) ] %>%
lapply(unlist) %>% lapply(min) %>% unlist()
subjectFeaturesList[[i]] <- cbind(gaitLength, gaitHeight, gaitDuration, expTime)
}
# Reading the segmented gait cycles for each participant:
# -------------------------------------------------------
profileFiles <- list.files(path = "../Data/matFiles/Profiles/", pattern = "*.mat", full.names = TRUE)
subjectProfilesList <- lapply(profileFiles, readMat) %>% # read all mat files from the path above
lapply("[[", "gait") # extract sublist titled "gait" from all the lists
names(subjectProfilesList) <- sapply(profileFiles, str_extract_all, "Subject[:digit:]*", simplify = TRUE)
# Reading task-related based on the tagging of the experimental videos:
# ---------------------------------------------------------------------
videoFiles <- list.files(path = "../Data/csvFiles/", pattern = "*.csv", full.names = TRUE)
videoDataList <- lapply(videoFiles, read.csv)
names(videoDataList) <- sapply(videoFiles, str_extract_all, "Subject[:digit:]*", simplify = TRUE)3 Data Preprocessing and Rational Subgroups
3.1 Data Preprocessing
In Section 2, we have read four different data sets and extracted the three gait features. In this subsection, we perform further preprocessing of the data. Specifically, we remove the first 10 minutes of experimental data from the videoDataList since these corresponded to the warm-up period of the experiment as detailed in Maman et al. (2017), Baghdadi et al. (2019), and Maman et al. (2020). WE NEED TO CONTINUE THE DISUCSSION OF WHAT WE HAVE DONE IN THIS SECTION. I STOPPED WORKING ON THIS UNTIL WE DISCUSS. Saeb, this is currently slightly different from your implementation as I am removing any task that spanned the warm-up period; I think this is a cleaner way of doing it and it will not change your preliminary results much
for (i in seq(length(rawSubjectsData))) {
videoDataList[[i]] %<>%
dplyr::select(Leaves, Enters) %>%
apply(MARGIN=2, FUN=na.omit) %>% # removing NA observations
as.data.frame %>%
mutate(Enters = Enters - 600, Leaves = Leaves - 600) %>% # subtracting 600 seconds for warm-up
filter(Enters > 0 & Leaves > 0) # removing observations/tasks, which contained any of the warm-up period
}3.2 Rational Subgroups
Given the nature of the experimental study, the subjects performed varying tasks that may affect our gait analysis; we have only examined data where the subjects are outside of the main laboratory room since this captures the material handling component of the examined task. An immediate consequence of using only “out-of-the-room” data is that the time differences between the retained segmented gait cycles is no longer almost constant since relatively large time gaps exist, which correspond to when the participant is in the room. For this reason, we have subgrouped the data. Consecutive “out-of-the-room” gait cycles are placed in one subgroup. Statistically speaking, our subgrouping strategy follows the rational subgrouping mantra in statistical surveillance. The purpose of rational subgrouping is two-fold:
- Minimize the variability within a subgroup such that the witin-group variability captures the natural variability in the gait cycles.
- Between group variability can be used to capture the effects of “assignable causes of variability”, which can include disease progression, fatigue, etc.
For more details on rational subgrouping, we refer the reader to Montgomery (2019) for more details. The code chunk below provides our approach to creating the rational subgroups for both the “profile” and “gait features” data sets.
numRows <- lengths(subjectProfilesList)/5
confidence <- videoDataList %>% lapply(function(x) {mean(x[,2] - x[,1])}) %>%
do.call(what=c) %>% `*` (0.05) # 5% confidence for each subgroup
subgroupList <- list()
for (id in seq(length(rawSubjectsData))) {
aM <- subjectProfilesList[[id]] %>% do.call(what=c) %>%
.[(numRows[[id]]+1):(2*numRows[id])] %>%
lapply(as.vector)
expT <- subjectProfilesList[[id]] %>% do.call(what=c) %>%
.[(4*numRows[[id]]+1):(5*numRows[id])] %>%
lapply(function(x) {x[1,1]}) %>% do.call(what=c)
###########################################################
# Match times of walking cycles, extracted from videos, against the experimental times
leaveIdx <- match.closest(videoDataList[[id]]$Leaves + confidence[id], expT)
enterIdx <- match.closest(videoDataList[[id]]$Enters - confidence[id], expT)
rmIdx <- c(3000, which(enterIdx - leaveIdx < 5))
leaveIdx <- leaveIdx[-rmIdx]
enterIdx <- enterIdx[-rmIdx]
########################################################## Create subgroups
aMList <- list()
LHDList <- list()
for (i in 1:length(leaveIdx)) {
aMList[[i]] <- aM[leaveIdx[i]:enterIdx[i]]
LHDList[[i]] <- subjectFeaturesList[[id]][leaveIdx[i]:enterIdx[i],]
}
######################################### Remove subgroups with length < 20
LHDList %<>% .[-{which(lengths(aMList) < 20) %>% c(3000)}]
aMList %<>% .[-{which(lengths(aMList) < 20) %>% c(3000)}]
######################################### Remove observations with profiles< 2nd percentile profile length; This is done to cut all profiles at this length (equal number of features)
cutLen <- aMList %>% do.call(what=c) %>% lengths() %>% quantile(probs=0.02)
for (i in 1:length(aMList)) {
rmIdx <- which(lengths(aMList[[i]]) < cutLen) %>% c(3000)
aMList[[i]] <- aMList[[i]][-rmIdx]
LHDList[[i]] <- LHDList[[i]][-rmIdx,]
}
######################################## Create a matrix with cutLen number of columns in each profile subgroup
aMList %<>% lapply(function(x) {do.call(rbind, x)}) %>%
lapply(function(x) {x[,1:cutLen]})
############################################ Lengths of subgroups 0&1: 10&20
subgroupList[[id]] <- list(Acc = list(X0=aMList[1:10], X1=aMList[11:30], X2=aMList[31:length(aMList)]),
LHD = list(X0=LHDList[1:10], X1=LHDList[11:30], X2=LHDList[31:length(LHDList)]))
}
names(subgroupList) <- c(paste0("subject", c(1:12, 14:15)))4 \(T^2\) Calculations
4.1 \(T^2\) Calculations on Gait Acceleration Profiles
Initialize.
T2Profiles <- list()
for (id in seq(length(subgroupList))) {
X0 <- subgroupList[[id]]$Acc$X0
############################# Phase0
############################# Remove single outliers in number of passes
pass <- 1
for (i in 1:pass) {
UCL0 <- BootstrapUCL0(x0=X0, alpha=0.05, iterations=200) # Find in Functions
######################### Find T2 for original X0 and remove outliers from subgroups
X0Pooled <- do.call(rbind, X0)
T2 <- X0Pooled %>% apply(FUN=mean, MARGIN=2) %>%
rep.row(n=nrow(X0Pooled)) %>% -(X0Pooled) %>%
apply(MARGIN=1, function(x) {x %*% ginv(cov(X0Pooled)) %*% x})
tmp <- mapply(function(x,y) {rep(y, x)}, x=sapply(X0, nrow), y=1:10) %>% do.call(what=c)
OCIdx <- T2 %>% split(f=tmp) %>% lapply(function(x) {c(which(x>UCL0), 3000)})
X0 %<>% mapply(FUN=function(x, y) {x[-y,]}, y=OCIdx) # Updated X0 after phase0
}
X0Pooled <- do.call(rbind, X0)
X0PooledMean <- colMeans(X0Pooled)
############################# Within batch covariance; computed from phase0 after outlier removal
S_w <- X0 %>% lapply(colMeans) %>%
mapply(FUN=rep.row, sapply(X0, nrow)) %>% # Mean of each subgroup
mapply(FUN=function(x,y) {x-y}, y=X0, SIMPLIFY=F) %>%
lapply(function(x) {t(x) %*% x}) %>%
Reduce(f=`+`) %>%
`/` (sum(lengths(X0)/ncol(X0[[1]])) - length(X0))
############################# Phase1&2
############################# Read phase 1&2 subgroups
X1 <- subgroupList[[id]]$Acc$X1
X2 <- subgroupList[[id]]$Acc$X2
######################### Compute T2 in phase1
X1Mean <- sapply(X1, colMeans) %>% t()
T2 <- rep.row(X0PooledMean, length(X1)) %>% -(X1Mean) %>%
apply(MARGIN=1, function(x) x %*% ginv(S_w) %*% x)
######################### Bootstrap UCL using phase 1 data
sampled <- sample(T2, 1000*length(T2), replace=T) %>% matrix(ncol=length(T2), byrow=T)
UCL <- apply(sampled, 1, quantile, probs=0.99) %>% mean()
######################### T2 for phase2 data
X12Mean <- sapply(c(X1, X2), colMeans) %>% t()
T2 <- rep.row(X0PooledMean, length(c(X1, X2))) %>% -(X12Mean) %>%
apply(MARGIN=1, function(x) x %*% ginv(S_w) %*% x)
T2Profiles[[id]] <- list(T2=T2, UCL=UCL)
}
names(T2Profiles) <- paste0("subject", c(1:12, 14:15))4.2 Plots of \(T^2\) on Gait Acceleration Profiles
Initialize.
for (id in seq(length(subgroupList))) {
T2 <- T2Profiles[[id]]$T2
UCL <- T2Profiles[[id]]$UCL
color <- ifelse(T2>UCL, "OC", "IC")
T2.df <- data.frame(subgroup=1:length(T2), T2=T2, col=color,
shape=c(rep("Phase 1", 20), rep("Phase 2", length(T2) - 20)))
T2Plot <- ggplot(data=T2.df) +
geom_point(aes(x=subgroup, y=T2, shape=shape, col=col), size=2.5) +
geom_line(aes(x=subgroup, y=T2)) +
geom_hline(yintercept=UCL) +
annotate("text", 2, UCL, vjust=-1, label="UCL", size=5) +
scale_shape_manual(values = c(6, 16)) +
scale_color_manual(values = c("black", "red")) +
ggtitle(paste0("Subject ", id)) +
theme_bw() +
theme(plot.title=element_text(size = 20, hjust=0.5, margin = margin(t=10, b=-30)),
axis.text=element_text(size=15),
axis.title=element_text(size=20),
legend.title=element_blank(), # This removes the legend title
legend.position=c(.1, .68),
legend.text=element_text(size=10),
legend.key = element_rect(fill="grey"))
cat('###',paste0("Subject", id), "{-}",'\n')
print(T2Plot)
cat('\n \n')
}Subject1
Subject2
Subject3
Subject4
Subject5
Subject6
Subject7
Subject8
Subject9
Subject10
Subject11
Subject12
Subject13
Subject14
4.3 RPE and Profiles \(T^2\) Charts
Initialize.
# Reading the RPE values from gait_hotelling GitHub Repo:
RPE <- read.csv("https://raw.githubusercontent.com/fmegahed/gait_hotelling/master/Data/csvFiles/RPE/RPE.csv") %>% .[,-14]
SubgroupTimeList <- subgroupList %>% lapply(function(x) {x$LHD[c(2,3)]}) %>%
lapply(function(x) {do.call(c, x)}) %>%
lapply(function(x) {sapply(x, function(a) {median(a[,4])})})
OCIdx <- lapply(T2Profiles, function(x) {which(x$T2>x$UCL)})
for (id in seq(length(SubgroupTimeList))) {
OCTime <- SubgroupTimeList[[id]] %>% .[OCIdx[[id]]]
RPETime <- list(which(RPE[,id+1]>=13) %>% min(), which(RPE[,id+1]>=15) %>% min()) %>%
sapply(function(x) {RPE[x,1]})
cat('###',paste0("Subject", id), "{-}",'\n')
plot(OCTime, rep(1, length(OCTime)), pch=16, xlim=c(0,11000), xlab="Experiment Time",
ylab="", main=paste0("Subject", id))
abline(v=RPETime, col=c("blue", "red"), lwd=2)
legend("topleft", legend=c(expression(RPE >= 13), expression(RPE >= 15)), col=c("blue", "red"),
lty=1, lwd=2)
cat('\n \n')
}Subject1
Subject2
Subject3
Subject4
Subject5
Subject6
Subject7
Subject8
Subject9
Subject10
Subject11
Subject12
Subject13
Subject14
4.4 \(T^2\) Calculations on Gait Features
Initialize.
T2Features <- list()
SWList <- list()
X0Updated <- list()
for (id in seq(length(subgroupList))) {
X0 <- subgroupList[[id]]$LHD$X0 %>% lapply(function(x) {x=x[,1:3]}) # Last column contains experiment time
############################# Phase0
############################# Remove single outliers in number of passes
pass <- 1
for (i in 1:pass) {
UCL0 <- BootstrapUCL0(x0=X0, alpha=0.05, iterations=200) # Find in Functions
######################### Find T2 for original X0 and remove outliers from subgroups
X0Pooled <- do.call(rbind, X0)
T2 <- rep.row(colMeans(X0Pooled), n=nrow(X0Pooled)) %>% -(X0Pooled) %>%
apply(MARGIN=1, function(x) {x %*% ginv(cov(X0Pooled)) %*% x})
tmp <- mapply(function(x,y) {rep(y, x)}, x=sapply(X0, nrow), y=1:10) %>% do.call(what=c)
OCIdx <- T2 %>% split(f=tmp) %>% lapply(function(x) {c(which(x>UCL0), 3000)})
X0 %<>% mapply(FUN=function(x, y) {x[-y,]}, y=OCIdx) # Updated X0 after phase0
}
X0Pooled <- do.call(rbind, X0)
X0PooledMean <- colMeans(X0Pooled)
############################# Within batch covariance; computed from phase0 after outlier removal
S_w <- X0 %>% lapply(colMeans) %>%
mapply(FUN=rep.row, sapply(X0, nrow)) %>% # Mean of each subgroup
mapply(FUN=function(x,y) {x-y}, y=X0, SIMPLIFY=F) %>%
lapply(function(x) {t(x) %*% x}) %>%
Reduce(f=`+`) %>%
`/` (sum(lengths(X0)/ncol(X0[[1]])) - length(X0))
############################# Phase1&2
############################# Read phase 1&2 subgroups
X1 <- subgroupList[[id]]$LHD$X1 %>% lapply(function(x) {x=x[,1:3]})
X2 <- subgroupList[[id]]$LHD$X2 %>% lapply(function(x) {x=x[,1:3]})
######################### Compute T2 in phase1
X1Mean <- sapply(X1, colMeans) %>% t()
T2 <- rep.row(X0PooledMean, length(X1)) %>% -(X1Mean) %>%
apply(MARGIN=1, function(x) x %*% ginv(S_w) %*% x)
######################### Bootstrap UCL using phase 1 data
sampled <- sample(T2, 1000*length(T2), replace=T) %>% matrix(ncol=length(T2), byrow=T)
UCL <- apply(sampled, 1, quantile, probs=0.99) %>% mean()
######################### T2 for phase2 data
X12Mean <- sapply(c(X1, X2), colMeans) %>% t()
T2 <- rep.row(X0PooledMean, length(c(X1, X2))) %>% -(X12Mean) %>%
apply(MARGIN=1, function(x) x %*% ginv(S_w) %*% x)
T2Features[[id]] <- list(T2=T2, UCL=UCL)
X0Updated[[id]] <- list(X0=X0)
SWList[[id]] <- list(S_w=S_w)
}
names(T2Features) <- paste0("subject", c(1:12, 14:15))
names(X0Updated) <- paste0("subject", c(1:12, 14:15))
names(SWList) <- paste0("subject", c(1:12, 14:15))4.5 Plots of \(T^2\) on Gait Features
Initialize.
for (id in seq(length(subgroupList))) {
T2 <- T2Features[[id]]$T2
UCL <- T2Features[[id]]$UCL
color <- ifelse(T2>UCL, "OC", "IC")
T2.df <- data.frame(subgroup=1:length(T2), T2=T2, col=color,
shape=c(rep("Phase 1", 20), rep("Phase 2", length(T2) - 20)))
T2Plot <- ggplot(data=T2.df) +
geom_point(aes(x=subgroup, y=T2, shape=shape, col=col), size=2.5) +
geom_line(aes(x=subgroup, y=T2)) +
geom_hline(yintercept=UCL) +
annotate("text", 2, UCL, vjust=-1, label="UCL", size=5) +
scale_shape_manual(values = c(6, 16)) +
scale_color_manual(values = c("black", "red")) +
ggtitle(paste0("Subject ", id)) +
theme_bw() +
theme(plot.title=element_text(size = 20, hjust=0.5, margin = margin(t=10, b=-30)),
axis.text=element_text(size=15),
axis.title=element_text(size=20),
legend.title=element_blank(), # This removes the legend title
legend.position=c(.1, .68),
legend.text=element_text(size=10),
legend.key = element_rect(fill="grey"))
cat('###',paste0("Subject", id), "{-}",'\n')
print(T2Plot)
cat('\n \n')
}Subject1
Subject2
Subject3
Subject4
Subject5
Subject6
Subject7
Subject8
Subject9
Subject10
Subject11
Subject12
Subject13
Subject14
4.6 RPE and Features \(T^2\) Charts
Initialize.
# Reading the RPE values from gait_hotelling GitHub Repo:
RPE <- read.csv("https://raw.githubusercontent.com/fmegahed/gait_hotelling/master/Data/csvFiles/RPE/RPE.csv") %>% .[,-14]
SubgroupTimeList <- subgroupList %>% lapply(function(x) {x$LHD[c(2,3)]}) %>%
lapply(function(x) {do.call(c, x)}) %>%
lapply(function(x) {sapply(x, function(a) {median(a[,4])})})
OCIdx <- lapply(T2Features, function(x) {which(x$T2>x$UCL)})
for (id in seq(length(SubgroupTimeList))) {
OCTime <- SubgroupTimeList[[id]] %>% .[OCIdx[[id]]]
RPETime <- list(which(RPE[,id+1]>=13) %>% min(), which(RPE[,id+1]>=15) %>% min()) %>%
sapply(function(x) {RPE[x,1]})
cat('###',paste0("Subject", id), "{-}",'\n')
plot(OCTime, rep(1, length(OCTime)), pch=16, xlim=c(0,11000), xlab="Experiment Time",
ylab="", main=paste0("Subject", id))
abline(v=RPETime, col=c("blue", "red"), lwd=2)
legend("topleft", legend=c(expression(RPE >= 13), expression(RPE >= 15)),
col=c("blue", "red"), lty=1, lwd=2)
cat('\n \n')
}Subject1
Subject2
Subject3
Subject4
Subject5
Subject6
Subject7
Subject8
Subject9
Subject10
Subject11
Subject12
Subject13
Subject14
5 Diagnosis of Gait Features
5.1 Decomposition of Features \(T^2\)
Initialize.
qList <- list()
for (id in seq(length(subgroupList))) {
############### Read inputs
X0 <- X0Updated[[id]]$X0
X1 <- subgroupList[[id]]$LHD$X1 %>% lapply(function(x) {x=x[,1:3]})
X2 <- subgroupList[[id]]$LHD$X2 %>% lapply(function(x) {x=x[,1:3]})
X0PooledMean <- do.call(rbind, X0) %>% colMeans()
X1Mean <- sapply(X1, colMeans) %>% t()
X12Mean <- sapply(c(X1, X2), colMeans) %>% t()
S_w <- SWList[[id]]$S_w
T2 <- T2Features[[id]]$T2
UCL <- T2Features[[id]]$UCL
t3 <- which(T2 > UCL)
p <- 3 # number of features
alpha <- 0.01
############## T2 Decomposition and UCL's
q <- init.decomp(p=p) # initialize the decomposition matrix
for(i in 1:nrow(q)) {# Control limits for each component of the decomposition
b <- subset(q[i,3:ncol(q)], q[i,3:ncol(q)] > 0)
T2 <- rep.row(X0PooledMean[b], length(X1)) %>% -(X1Mean[,b]) %>%
apply(MARGIN=1, function(x) x %*% solve(S_w[b,b]) %*% x)
sampled <- matrix(sample(T2, 1000 * length(T2), replace=T), ncol=20, byrow=T)
q[i,2] <- apply(sampled, 1, quantile, probs=1-alpha) %>% mean()
}
############## Decompose T2 for every out of control data point
tmpList <- list()
for(ii in 1:length(t3)){
for(i in 1:nrow(q)) {
b <- subset(q[i,3:ncol(q)], q[i,3:ncol(q)] > 0)
q[i,1] <- t(X0PooledMean[b]-X12Mean[t3[ii],][b]) %*%
solve(S_w[b,b]) %*%
(X0PooledMean[b]-X12Mean[t3[ii],][b])
}
############ Store decomposed elements
tmpList[[ii]] <- q
}
names(tmpList) <- paste0("subgroup ", t3)
qList[[id]] <- tmpList
}
names(qList) <- paste0("subject", c(1:12, 14:15))5.2 2D Ellipse plots
Initialize.
################ Select subject number and features
id <- 1 # Select the subject number
f <- c(2,3) # Select 2 features by interpreting the T2 decomposition
################ Read inputs
X0 <- X0Updated[[id]]$X0
X1 <- subgroupList[[id]]$LHD$X1 %>% lapply(function(x) {x=x[,1:3]})
X2 <- subgroupList[[id]]$LHD$X2 %>% lapply(function(x) {x=x[,1:3]})
X0Pooled <- do.call(rbind, X0)
X0PooledMean <- do.call(rbind, X0) %>% colMeans()
X1Pooled <- do.call(rbind, X1)
S_w <- SWList[[id]]$S_w
q <- qList[[id]][[1]] # Read one of the q's as its 2nd column (q[,2]) are the control limits and thus the same across all subgroups in one subject
T2 <- T2Features[[id]]$T2
UCL <- T2Features[[id]]$UCL
################ Univariate CL's and 2D plot limits
XLim <- c(X0PooledMean[f[1]]-sqrt(S_w[f[1],f[1]]*q[f[1],2]),
X0PooledMean[f[1]]+sqrt(S_w[f[1],f[1]]*q[f[1],2]))
YLim <- c(X0PooledMean[f[2]]-sqrt(S_w[f[2],f[2]]*q[f[2],2]),
X0PooledMean[f[2]]+sqrt(S_w[f[2],f[2]]*q[f[2],2]))
xlim <- c(min(XLim) - 1.2*diff(XLim), max(XLim) + 1.2*diff(XLim))
ylim <- c(min(YLim) - 1.2*diff(YLim), max(YLim) + 1.2*diff(YLim))
################ Create grid for the ellipse contour plot
ngrid <- 100
grid <- cbind(seq(xlim[1], xlim[2], length = ngrid),
seq(ylim[1], ylim[2], length = ngrid))
################ T2 on the grid points
T2_grid <- apply(expand.grid(grid[,1], grid[,2]), 1,
function(x) {(x- X0PooledMean[f]) %*% solve(S_w[f,f]) %*% (x- X0PooledMean[f])})
T2_grid <- matrix(T2_grid, ngrid, ngrid)
######################################################## Plot
############## Plot the contour
# old.par <- par()
par(mar=c(5, 5, 4, 2) + 0.1)
plot(1 ,1, type = "o", xlim = xlim, ylim = ylim, main="Ellipse Chart", cex.main=1.5,
xlab=paste0(colnames(X0[[1]])[f[1]], " (m)"),
ylab=paste0(ylab=colnames(X0[[1]])[f[2]], " (s)"),
cex.axis=1.5, cex.lab=1.5)
legend("topleft", legend=c("phase 0 subgroup means", "phase 1 subgroup means",
"X control limits", "Y control limits"),
col=c("green", "blue", "black", "brown"), pch=c(3,3,NA,NA), lty=c(NA,NA,2,2),
cex=1, pt.cex=2, lwd=2, bg="grey")
################
idx <- apply(q[,3:5], MARGIN=1, function(r) all(c(f,0) %in% r)) %>% which()
contour(grid[,1], grid[,2], T2_grid, levels=q[idx,2], drawlabels=T, add=T)
################ Plot phase 0 data
center0 <- X0 %>% lapply(function(x) {x[,f]}) %>% sapply(colMeans) %>% t()
points(center0[,1], center0[,2], pch=3, cex=2, col="green", lwd=2)
# points(X0Pooled[,f], pch=1, cex=0.5, col="blue")
################ Plot phase 1 data
center1 <- X1 %>% lapply(function(x) {x[,f]}) %>% sapply(colMeans) %>% t()
points(center1[,1], center1[,2], pch=3, cex=2, col="blue", lwd=2)
# points(X1Pooled[,f], pch=1, cex=0.5, col="blue")
################ plot phase 1&2 data by subgroup number
center12 <- c(X1,X2) %>% lapply(colMeans) %>% sapply(function(x) {x[f]}) %>% t()
labels <- T2 %>% t() %>% as.data.frame() %>% names() %>% gsub(pattern='V', replacement='')
text(center12, labels=labels, cex = 1, col=ifelse(T2>UCL, "red", "black"))
################ Plot univariate CL's
abline(v=XLim, lty="dashed", col="black", lwd=2)
abline(h=YLim, lty="dashed", col="brown", lwd=2)
References
Baghdadi, A., Cavuoto, L. A., Jones-Farmer, A., Rigdon, S. E., Esfahani, E. T., & Megahed, F. M. (2019). Monitoring worker fatigue using wearable devices: A case study to detect changes in gait parameters. Journal of Quality Technology, Online First. https://doi.org/10.1080/00224065.2019.1640097
Maman, Z. S., Chen, Y.-J., Baghdadi, A., Lombardo, S., Cavuoto, L. A., & Megahed, F. M. (2020). A data analytic framework for physical fatigue management using wearable sensors. Expert Systems with Applications, 113405.
Maman, Z. S., Yazdi, M. A. A., Cavuoto, L. A., & Megahed, F. M. (2017). A data-driven approach to modeling physical fatigue in the workplace using wearable sensors. Applied Ergonomics, 65, 515–529.
Montgomery, D. C. (2019). Introduction to statistical quality control. Wiley.
Email: saebraga@buffalo.edu | Website: LinkedIn↩︎
Email: jiyeonk@buffalo.edu | Phone: +1-716-645-6063 | Website: University at Buffalo Official↩︎
Email: loracavu@buffalo.edu | Phone: +1-716-645-4696 | Website: University at Buffalo Official↩︎
Email: fmegahed@miamioh.edu | Phone: +1-513-529-4185 | Website: Miami University Official↩︎
Email: farmerl2@miamioh.edu | Phone: +1-513-529-4823 | Website: Miami University Official↩︎